Self-consistent Single-band Approximation for Interacting Boson Systems 
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Traditionally, the single-band approximation for interacting many-body systems is done with 
pre-determined single-particle Wannier functions, ignoring the dependence of the Wannier function 
on interaction. We show that the single-band approximation has to be done self-consistently to 
properly account the interaction effect on the Wannier functions. This self-consistent single-band 
approximation leads to a nonlinear equation for Wannier functions, which can be recast into a set 
of nonlinear equations for Bloch functions. These equations are simplified for two special cases, 
the superfluid regime and deep in the Mott insulator regime. A simple example with double-well 
potential is used to illustrate our results. 
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Interacting many-body systems are notoriously hard 
to solve. One usual tactic is the single-band approxima- 
tion, where the Wannier functions [l[ for the lowest Bloch 
band are used to reduce the system to a lattice model, 
such as Hubbard model @, OJ. Traditionally, the Wan- 
nier functions used in the approximation are obtained ei- 
ther directly from the single-particle Bloch waves or with 
some variational approaches 0, The interaction effect 
on the Wannier function is completely neglected in all 
these traditional methods. This neglect is at odds with 
the fact that the interaction evidently affects the shape 
of the Wannier functions. 

Current strong interest in ultracold atoms in optical 
lattices has put this problem into spotlight Q. A recent 
experiment with ultracold atoms clearly demonstrated 
that the on-site interaction, therefore, the Wannier func- 
tions, depends on the occupation number 0- It has also 
seen increased theoretical efforts to address this problem 
@; H [3 U 1 13, 3 ■ However, all the efforts have certain 
drawbacks from a general perspective. In the variational 
approach adopted by Li et al. 0], there is arbitrariness 
in the choice of the trial Wannier function. The mean- 
field method used by Liang et al. can apply only in the 
superfluid regime and can not be justified to apply for 
the more interesting superfluid-Mott insulator transition 
regime (loj. The ab initio calculation in Ref. [ll| is hard 
to be scaled up for systems with large number of parti- 
cles. 

We address the interaction effect on Wannier function 
by re-examining the single-band approximation. We find 
that the interaction effect can be taken into account au- 
tomatically once the single-band approximation is done 
self-consistently. The self-consistent single-band approx- 
imation is achieved by minimizing the ground state en- 
ergy of the system by varying Wannier function. The 
minimization leads to a nonlinear equation for the Wan- 
nier function, which depends on the ground state of the 
system. Therefore, one has to solve the nonlinear equa- 
tion self-consistently with the ground state of the re- 
sulted lattice model to find the interaction-dependent 



Wannier functions. Our results provide a general varia- 
tional framework for transforming a periodic system into 
a lattice model self-consistently and can be generalized 
to systems where multiple bands are needed. 

One can reformulate the nonlinear equation for Wan- 
nier functions in terms of Bloch functions, and obtain a 
set of nonlinear equations for the Bloch functions. Sim- 
plified forms for these nonlinear equations are obtained 
for two special cases, the superfluid regime and deep in 
the Mott insulator regime. Our results are illustrated 
with a system of double-well potential, where some gen- 
eral properties of the nearest neighbor tunneling param- 
eter J and the on-site interaction U are revealed. 

Although our approach can be applied to other sys- 
tems, e.g., fermionic systems, we here focus on the sys- 
tem of ultracold bosons in optical lattices [H, 0, EH- For 
bosons of mass m, the Hamiltonian is 



H 



drip* (r) 



2m 



V(r) ^(r) 



drdr' 



^(v)^(r')U(\r~v'\)^(v') , (1) 



where v(r) is a periodic potential and u(|r|) is the inter- 
action between two atoms. If one is interested only in the 
low temperature properties of the system, one can use the 
single-band approximation and expand the bosonic field 
operator as 



^(r)=^a i Wj(r), 



(2) 



where Wj(r) 
site j and a. 



= W(r — Yj) is the Wannier function at 
is the associated annihilation operator. 
In this case, the trial ground state \G t ) is given by 
\G t ) = F(a])|vaccum) , where the functional F is to be 
determined. Usually, the Wannier function in Eq.([2]) is 
pre-determined. Here the Wannier function is not known 
a priori except that it is expected to resemble the single- 
particle Wannier function. We look for the Wannier 
functions that minimize the system's single-band ground 



2 



state \Gt 



E G = (G t \H\G t ) = (G t \H bh \G t ), 



where Hhh is the usual Bose-Hubbard Hamiltonian 



333A 



H bh = -y] Jhh^j^h + ^2 u ji 



The parameters are given by 

J ili2 = - J dvWl{v)H W h {v) 



U, 



drdr' 



W£(r)W£(r^ 



xU(\r-r'\)W j3 (r)W u (r') 



(3) 



(4) 



(5) 



(6) 



where Hq = — |^-V 2 + V(v) . Usually further approxima- 
tion is made so that only two terms, the nearest neigh- 
bor tunneling and the on-site interaction, are kept in the 
Bose-Hubbard model. For the purpose of deriving gen- 
eral formalism, this approximation is not needed so far. 

We achieve the minimization of the ground state en- 
ergy Eq by varying the Wannier function under the or- 
thonormal constraints hj = / drW* (r)W(r — rj) = 6o,j ■ 
According to the Feynman-Hellman theorem, we have 



6{E G - Ej MA) 
SW*(r) 



N 5h i 
SW*(r) 



= 0, 



(7) 



where /i's are Lagrangian multipliers, we obtain a non- 
linear equation for the Wannier functions 



3 



31,32 

3334 



E<4 



3132 



dr' 



W*(r' 



r 3l) 



W(r' + r h - r j3 )U(\r' - r|) W(r + r h - r j4 ) 



(8) 



It is clear that the above equation depends on the ground 
state \Gt) of the Bose-Hubbard model in Eq.(|l]), which 
can be found with many-body methods such as direct 
diagonalization, Gutzwiller proiection(l6j. density ma- 
trix renormalization group (DMRG)[l7l, or time-evolving 
block decimation (TEBD)[l8|. At the same time, we 
need the Wannier functions to compute J's and t/'s for 
the Bose-Hubbard model. Therefore, the Bose-Hubbard 



model in Eq.Q has to be solved self-consistently with 
the above nonlinear equation. We call this self-consistent 
single-band approximation. Due to the complexity of the 
equations, one can apply further approximations. For ex- 
ample, one can opt to completely ignore the interaction 
while solving Eq.©. This is just what people have tradi- 
tionally done with the single-band approximation. One 
can also choose to keep only the nearest neighbor tunnel- 
ing J and the on-site interaction U in the Bose-Hubbard 
Hamiltonian (H}. However, the off-site terms in the last 
summation in Eq. |S]) can not be dropped simultaneously. 

A periodic system can be described alternatively with 
Bloch functions. If we place the system in a box of N 
lattice sites, the Wannier functions are related to Bloch 
functions as 



W(r - r„) 



THE- 



— ik-r 



*k(r). 



(9) 



where is a Bloch function with Bloch wave number 
k and is normalized to one. For Bloch functions, the 
nonlinear equation JSJ) becomes 

<kik 2 



x / dr' 



* ki (rO* k3 (r'M|r'-r|) * k4 (r) 



(10) 



where <kik2k3k4> stands for summation with the con- 
straint ki+k 2 = k 3 +k 4 +K, 6 k = ^= J2 n d n e~ lk rn , and 

= j? J2 n A*n e_lk r ™ ■ K is a reciprocal lattice. While 
we can not split Eq. ((8]) into a set of equations for Wan- 
nier functions at different sites, we are allowed to split 
the above equation for different Bloch wave numbers k 
and obtain 

<k x k 



k 3 k 4 > 

* k3 (r') x U(\r'-v\) * k4 (r) 



£k*k(r) - ^o*k(r)+ ]T P kl kk 3 k 4 / rfr'[* ki (r') 

(11) 



where P klkk3k4 = (& kl & k &k 3 &k 4 )/(& k & k) and £ k = 
f k /(S k 6 k ). In the following discussion, for simplicity, we 
shall use dilute atomic gases, where M(|r|) = goS(r), [l9[ 
to discuss two special cases. 

We consider first the superffuid regime. With the Bo- 
goliubov mean- field theory [l9j |. we have 



bh 



k^O 

X)M)^k(S k Sl k + SkS. 

k#0 



(12) 



where jVq is the number of atoms in the state ipo, £ k 
J dnfcHofa, and U k = (g /2) J dr|V k | 2 |^o| 2 - Following 
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the standard procedure [19j , we find that 

-Pk,kk,k, = Ufc 1 (^ki,k 4 ^k,k 3 + ^ki,k 3 ^k,k 4 ) 

-5-k 1 ,k5-k 3 ,k 4 - 2Jki,0<5k,oA/o , (13) 



kikk 3 k 4 

U* k Uk 3 Vk 3 



where uq = vq 
and w k -^Q = 1 - 



,,2 



k^O 



[(e k + 4Wk)/£ k -l]/2 , 
V(£k + W ;7 k ) 2 -4A^^. 
This leads to a set of simplified nonlinear equations for 
Bloch functions 



i>o*o = i?o *o + ffaA/o | *o 1 2 *o 



5o5Z u kl*k| 2 *o 

k^O 



(14) 



and for 
^k*k = 



#o*k 



9o 



25oA^ |*o| 2 *k 
V \Uk>v k t— + 2v\, 



k'^0 



l*k'r*i, 



(15) 



Since and i>fc themselves depend on ^ k , the above two 
equations have to be solved self-consistently. Note that 
in the above derivation we have assumed that the lattice 
potential is symmetric, V(r) = V{— r), so that V>-k = V'k- 
We have also ignored the scattering processes with non- 
zero K. 

In contrast to superfluid regime, deep in the Mott- 
insulator regime, we have 



— (n + no)Sj 1 1„ <5 10 i,(5-i' 



31 J2 U J2 ,33 u 33 ,34 



where no is the averaged number of bosons per site. In 
this case, Eq.© is simplified and has the form 



£w(r) = H W(r) 



+go(n -l)\W(r)\ 2 W(r) 



(18) 



The off-site terms on the right hand side is usually very 
small and can be ignored. With this in mind, we imme- 
diately have one observations. With one atom per site, 
i.e., hq = 1, the Wannier function is approximately the 
single-particle ground state wave function of each indi- 
vidual well of the lattice; there is no interaction effect on 
the Wannier function. 

We now use one dimensional double- well potential with 
periodic boundary condition to illustrate our theory. The 
two Wannier functions for the left well and the right well 
are related to the ground state and the first excited state 
for the double-well potential as follows 
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FIG. 1: Wannier functions for different depths of the double- 
well potential. For clarity, the curves for v = 9.0 and v = 15.0 
are shifted up by 0.5 and 1.0, respectively, c = 2.5. The unit 
of v is the recoil energy E r = Ti 2 k\/2m. 



where ^0 and \&i are chosen so that they are both pos- 
itive in the left well. These two Wannier functions sat- 
isfy the nonlinear equation 1(5)1 . The corresponding Bose- 
Hubbard model is 

H 2 = -J + 2(No-l)U 3 {aja r + alai) + 

+ U2(a\,a\.aiai + Aa[aia\.a r + a\a\a r a r ) 

+ U(ajajaiai + a\.a\.a r a r ) , (20) 



where U 2 = Uiirr, U 3 = Ui rrri and U = Uuu- 

We choose the double-well potential as a part of the 

(16) one dimensional optical lattice created experimentally 
in Ref . [20l| . So, the double well potential is given by 
Vix) — Vosin 2 (fcix), where is the wave number of 

(17) the laser that creates the potential. Due to the lat- 
eral confinement, the interaction strength g is given 
by go = (4:nh 2 a s /rn)(muj±/(27rh)) — 2huj±a s , where 
a s is the s-wave scattering length and u± the perpen- 
dicular confinement frequency. We consider the case 

Tj)\ 2 W(r) of no = 1(Nq = 2). In our numerical calculations, 

the Hamiltonian in Eq. (j20|l is diagonalized directly and 
Eq.® is solved with the nonlinear equation solver in 
MATLAB. And we use c = imgo/(h ki,) as the dimen- 
sionless interaction parameter. 

The Wannier functions found numerically are plotted 
in FigJTJ As expected, the Wannier function becomes 
more and more localized as the depth of the well in- 
creases. It is worthwhile to note that the Wannier func- 
tions for shallow wells (e.g., v = 3.0) have nodes, which 
shows that nodeless ground state does not necessarily 
imply a nodeless Wannier function. 

The numerical results for the dependence of J and U 
on the well depth and interaction strength are shown in 
terms of the ratios J/ Jo and U/Uq in Fig|2] J and Uq 
are the tunneling parameter and on-site interaction ob- 
tained with single-particle Wannier function. Certain in- 
teresting behaviors of J/ Jo and U /Uq are revealed. For a 

(19) 

fixed interaction strength, both J/ Jo and U/Uq approach 
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FIG. 2: The tunneling parameter J and the on-site interaction 
Ui as functions of the double- well depth. For (a) and (c), 
c = 2.5; for (b) and (d), v = 5.0. The energy unit is the recoil 
energy E r — h 2 k 2 L /2m. 



one as the lattice gets stronger. If the lattice strength 
is fixed and the interaction strength changes, J/ Jo and 
U/Uo reach their respective extremum values at a critical 
interaction strength. 

We argue that the behaviors of J/ Jo and U/Uo re- 
vealed in Figgis not limited to the double- well systems, 
and should hold generally for ultracold bosons in optical 
lattices. As already discussed, deep in the Mott insulator 
regime, the solution of Eq. (fT5]) is single-particle Wannier 
function for n = 1. This implies that both of the ratios 
J I Jo and U/Uq should approach one as the optical lat- 
tice gets stronger with fixed interaction strength. This is 
exactly what is seen in Fig. [2a, c). If one fixes the lattice 
strength and increases the repulsive interaction, then one 
has single-particle Wannier function at both the begin- 
ning (c = 0) and the end (c» 1) of this change. This 
indicates that J/ Jo should reach its maximum point and 
U/Uo arrive at its minimum point during this process as 
shown in Fig. [5Jb,d). If one simulates the superfiuid- 
Mott insulator transition with the self-consistent single 
band approximation proposed here, the behaviors of J / Jo 
and U/Uq m Fig- EIb,d) can be used to determine the 
transition point between superfluid and Mott insulator. 

Note that the interaction effect on J and U is much 
smaller compared to the mean- field results in Ref . [Tol| . 
This is likely due to the finite size of our example sys- 
tem, where the Wannier function is confined and can not 
extend to infinity. More extensive numerical studies are 
needed to clarify this. If future numerical computation 
indeed indicates that the interaction effect on J and U is 
as large as indicated in the mean-field results [loj] . then 
the behaviors of U/Uo m Fig. [2c, d) may be used as an 
experimental means to detect the superfluid-Mott insula- 
tor transition as U can now be measured experimentally 
0| and the interaction can be adjusted with the Feshbach 
resonance 21 1 . 



Currently in typical experiments, ultracold atoms are 
also trapped by a harmonic potential 14 , [l5| . Conse- 
quently, the wells are not identical to each other. Also a 
random potential can be added to make the wells non- 
identical. Nevertheless, single-band approximation can 
still be applied as long as the difference between the 
site energies is smaller than the energy gaps between the 
ground state and the first excited state in the wells. For 
simplicity, we consider a one dimensional potential of N 
wells, which are not identical. In this case, the Wannier 
functions for the lowest "band" can be defined as 



1 



v k=0 



N-l 

^ — -V - 2kjK 



j = 0,l,r--,N-l, (21) 



where V>fc's are the lowest N eigenstates. Our variational 
approach can be easily adopted to this case with just one 
modification. Since the Wannier functions at different 
sites have different shapes, the constraint is now h n , m = 
f dxW*(x)W m (x) = 5 n>m . As a result, one obtains a set 
of nonlinear equations for the Wannier functions. Since 
everything is straightforward, we shall not write out the 
equations here. 
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